function F = solve_eqm_cf1_inner(param, const, integ, eqm)
%% Counterfactual 1: Baseline

    % normalization
    P_s     = param.P_s; 
    e       = param.e; 

    % read in parameters
    kappa   = param.kappa;
    alpha_m = param.alpha_m; 
    alpha_s = param.alpha_s;
    beta_v  = param.beta_v; 
    beta_m  = param.beta_m;    
    sigma   = param.sigma; 
    %sigma_s = param.sigma_s;
    f_mH    = param.f_mH;
    f_mF    = param.f_mF;
    f_vH    = param.f_vH;
    f_vF    = param.f_vF;
    mu_E    = param.mu_E;
    sigma_E = param.sigma_E;  
    mu_I    = param.mu_I;
    sigma_I = param.sigma_I;
    w       = param.w;
    Q       = param.Q;
    Nf      = param.Nf;
    
    % read in constants
    gamma   = const.gamma;
    k0      = const.k0;
    C_m0    = const.C_m0;
    C_v0    = const.C_v0;
    C_x00    = const.C_x00;   
    phi_v   = const.phi_v;
    phi_y   = const.phi_y;
    phi_s   = const.phi_s;
    D_F     = const.D_F;
    c_F     = const.c_F;
    lnzQ    = const.lnzQ;
        
    % read in all the integral from previous round
    E_zm00 = integ.E_zm00;
    E_zm01 = integ.E_zm01;
    E_zm10 = integ.E_zm10;
    E_zm11 = integ.E_zm11;
    E_zv00 = integ.E_zv00;
    E_zv01 = integ.E_zv01;
    E_zv10 = integ.E_zv10;
    E_zv11 = integ.E_zv11;
    E_zx00 = integ.E_zx00;
    E_zx01 = integ.E_zx01;
    E_zx10 = integ.E_zx10;
    E_zx11 = integ.E_zx11;
        
    % adjust for level in case integrals too large
%     rescale = max(1,10^floor(log10(max(E_zx11))-1));
%     E_zx00 = E_zx00/rescale;
%     E_zx01 = E_zx01/rescale;
%     E_zx10 = E_zx10/rescale;
%     E_zx11 = E_zx11/rescale;
%     E_zv00 = E_zv00/(rescale^(1/beta_v));
%     E_zv01 = E_zv01/(rescale^(1/beta_v));
%     E_zv10 = E_zv10/(rescale^(1/beta_v));
%     E_zv11 = E_zv11/(rescale^(1/beta_v));
%     E_zm00 = E_zm00/(rescale^(1/beta_m));        
%     E_zm01 = E_zm01/(rescale^(1/beta_m));        
%     E_zm10 = E_zm10/(rescale^(1/beta_m));        
%     E_zm11 = E_zm11/(rescale^(1/beta_m));        

    
%% Solve the equilibrium    

    % Setting
    tol = 1e-12;
    maxIter = 1000;

    % Initial guess
    D_Hnew = eqm.D_H; 
    c_Hnew = eqm.c_H;

    % Iteration
    err = 1;
    iter = 0;
    while err >= tol && iter < maxIter

        % new guess
        D_H = D_Hnew;
        c_H = c_Hnew;

        % D(q,E)
        D0 = D_H;
        rv_ads = (D_H/f_vH).^(1/(beta_v-1))./( (D_H/f_vH).^(1/(beta_v-1))+(e^sigma*D_F/f_vF).^(1/(beta_v-1)) );
        D1 = rv_ads.*D_H + (1-rv_ads)*e^sigma.*D_F;        
        
        % c(q,I)
        c0 = c_H;
        rm_ads = (c_H.^(1-sigma)/f_mH).^(1/(beta_m-1))./( (c_H.^(1-sigma)/f_mH).^(1/(beta_m-1)) + ((e*c_F).^(1-sigma)/f_mF).^(1/(beta_m-1)) );
        c1 = ( rm_ads.*c_H.^(1-sigma) + (1-rm_ads).*(e*c_F).^(1-sigma) ).^(1/(1-sigma));

        % f_v(q,E), f_m(q,I)
        f_v1 = rv_ads.^beta_v*f_vH + (1-rv_ads).^(beta_v)*f_vF;
        f_m1 = rm_ads.^beta_m*f_mH + (1-rm_ads).^(beta_m)*f_mF;

        % C_v(q,E), C_m(q,I) 
        C_v1 = f_v1.^(-1/beta_v).*(1./(sigma*w)).^(1/beta_v);
        C_m1 = f_m1.^(-1/beta_m).*(alpha_m./(sigma*w)).^(1/beta_m);

        % C_x(q,E,I)
        tempterm = (param.sigma/(param.sigma-1)*param.w.^(1-param.alpha_m-param.alpha_s)).^(1-param.sigma);
        C_x11 = (tempterm.*C_m1.^param.alpha_m.*C_v1).^const.gamma;
        C_x10 = (tempterm.*C_m0.^param.alpha_m.*C_v1).^const.gamma;
        C_x01 = (tempterm.*C_m1.^param.alpha_m.*C_v0).^const.gamma;         

        % Pi(q,E,I)
        Pi11 = D1.^gamma.*(c1.^alpha_m*P_s^alpha_s).^((1-sigma)*gamma).*C_x11;
        Pi10 = D1.^gamma.*(c0.^alpha_m*P_s^alpha_s).^((1-sigma)*gamma).*C_x10;
        Pi01 = D0.^gamma.*(c1.^alpha_m*P_s^alpha_s).^((1-sigma)*gamma).*C_x01;     
        Pi00 = D0.^gamma.*(c0.^alpha_m*P_s^alpha_s).^((1-sigma)*gamma).*C_x00;
      
        % M(q)
        M =         C_m0.*(Pi00.^(1/beta_m).*E_zm00 + Pi10.^(1/beta_m).*E_zm10) ...
          + rm_ads.*C_m1.*(Pi01.^(1/beta_m).*E_zm01 + Pi11.^(1/beta_m).*E_zm11); 
        
        % V(q)
        Vbar =         C_v0.*(Pi00.^(1/beta_v).*E_zv00 + Pi01.^(1/beta_v).*E_zv01)...
             + rv_ads.*C_v1.*(Pi10.^(1/beta_v).*E_zv10 + Pi11.^(1/beta_v).*E_zv11);      
        V = sum(phi_v.*repmat(Vbar,Q,1),2)';

        % theta_v(q) and theta_m(q)
        xi = M./V;
        theta_v = 1-exp(-kappa*xi);
        theta_m = (1-exp(-kappa*xi))./xi;
        theta_m(isnan(theta_m))=kappa;       
        theta_m(theta_m==0)=kappa;
        
        % X(q,E,I)
        X00 = Pi00.*E_zx00;
        X01 = Pi01.*E_zx01;
        X10 = Pi10.*E_zx10;
        X11 = Pi11.*E_zx11;
       
        % X_m(q)
        rm_spend = rm_ads.*c_H.^(1-sigma)./(c1.^(1-sigma));
        X_m = alpha_m*(sigma-1)/sigma*((X10+X00)+rm_spend.*(X11+X01));

        % P(q)^(1-sigma)
        Psig = (X00+X01)./D0 + rv_ads.*(X10+X11)./D1;
        
        % D_m(q)
        TempD = theta_v./M.*c_H.^(sigma-1).*X_m;
        TempD(isnan(TempD))=0;
        D_m = sum(phi_y.*phi_v.*repmat((TempD)',1,Q),1);

        % manufacturing trade blance: exports - imports
        rv_rev = rv_ads.*D_H./D1;
        B1 = sum((1-rv_rev).*(X11+X10)) - alpha_m*(sigma-1)/sigma*sum((1-rm_spend).*(X11+X01));

        % X_s: initial equilibrium level
        X_s = eqm.X_s;

        % D_s(q)
        D_s = phi_s.*X_s/sum(phi_s.*Psig);

        % D(q) (new guess)
        D_Hnew = D_m+D_s;

        % c(q) (new guess)
        c_Hnew = (theta_m./V.*sum(phi_y.*phi_v.*repmat(Psig,Q,1),2)').^(1/(1-sigma));   
        
        % update
        D_Hnew = D_H + 0.2*(D_Hnew - D_H);
        c_Hnew = c_H + 0.2*(c_Hnew - c_H);
        iter = iter+1;      

        % check convergence
        err1 = max(abs(D_Hnew-D_H)./D_H);
        err2 = max(abs(c_Hnew-c_H)./c_H);
        err = err1 + err2;
        
    end
    
    % Warning
    if iter == maxIter
        msg = 'User Warning: Inner Loop Maximum Iterations Reached.';
        warning(msg)
    end 

    % Adjust back
%     Pi00 = Pi00/rescale;
%     Pi01 = Pi01/rescale;
%     Pi10 = Pi10/rescale;
%     Pi11 = Pi11/rescale;
    

%% Additional

    % lnZ(omega,q)
    lnZQ = k0*lnzQ-log(sigma*gamma);

    % Import Decision
    ImportCutoffQ0 = (lnZQ+log(kron(ones(Nf,1),Pi01)-kron(ones(Nf,1),Pi00))-mu_I)/sigma_I;
    ImportCutoffQ1 = (lnZQ+log(kron(ones(Nf,1),Pi11)-kron(ones(Nf,1),Pi10))-mu_I)/sigma_I;
    ImportProbQ0 = normcdf(ImportCutoffQ0);   
    ImportProbQ1 = normcdf(ImportCutoffQ1);   

    % Export Decision
    EPiQ0 = exp(lnZQ).*(kron(ones(Nf,1),Pi00) + kron(ones(Nf,1),Pi01-Pi00).*ImportProbQ0) ...
        - exp(mu_I+1/2*sigma_I^2).*normcdf(ImportCutoffQ0-sigma_I);
    EPiQ1 = exp(lnZQ).*(kron(ones(Nf,1),Pi10) + kron(ones(Nf,1),Pi11-Pi10).*ImportProbQ1) ...
        - exp(mu_I+1/2*sigma_I^2).*normcdf(ImportCutoffQ1-sigma_I);
    NetEPiQ = EPiQ1-EPiQ0;
    %min(min(NetEPiQ))
    if min(min(NetEPiQ))<0
        disp('Inner loop adjustment: EPiQ1-EPiQ0 < 0')
        NetEPiQ(NetEPiQ<0) = 0;
    end  
    ExportCutoffQ = (log(NetEPiQ)-mu_E)./sigma_E;
    ExportProbQ = normcdf(ExportCutoffQ);   

    % Save mbar (in service price)
    %P_sF = P_s/e*(B1/(X_s+B1))^(1/(1-sigma_s));    
    %mbar = sum(V)/sum(phi_s.*Psig)*(P_s^(1-sigma_s)-(e*P_sF)^(1-sigma_s))^((1-sigma)/(1-sigma_s));
    %f = X/(sigma*gamma*P_s);
    

%% Return

    eqm_cf.Pi00 = Pi00;
    eqm_cf.Pi01 = Pi01;
    eqm_cf.Pi10 = Pi10;
    eqm_cf.Pi11 = Pi11;
    eqm_cf.D_H = D_H;
    eqm_cf.D_F = D_F;
    eqm_cf.D0 = D0;
    eqm_cf.D1 = D1;
    eqm_cf.D_s = D_s;
    eqm_cf.D_m = D_m;
    eqm_cf.c_H = c_H;
    eqm_cf.c_F = c_F;
    eqm_cf.c0 = c0;
    eqm_cf.c1 = c1;
    eqm_cf.Psig = Psig;
    eqm_cf.P = Psig.^(1/(1-sigma));
    eqm_cf.X_m = X_m;
    eqm_cf.X00 = X00;
    eqm_cf.X01 = X01;
    eqm_cf.X10 = X10;
    eqm_cf.X11 = X11;
    eqm_cf.V = V;
    eqm_cf.Vbar = Vbar;
    eqm_cf.M = M;
    eqm_cf.xi = xi;
    eqm_cf.theta_v = theta_v;
    eqm_cf.theta_m = theta_m;  
    eqm_cf.rv_ads = rv_ads;
    eqm_cf.rv_rev = rv_rev;
    eqm_cf.rm_ads = rm_ads;
    eqm_cf.rm_spend = rm_spend;
    eqm_cf.C_v1 = C_v1;    
    eqm_cf.C_m1 = C_m1;    
    eqm_cf.C_x11 = C_x11;    
    eqm_cf.C_x10 = C_x10;    
    eqm_cf.C_x01 = C_x01;    
    eqm_cf.B1 = B1;
    eqm_cf.P_s = P_s;
    eqm_cf.e = e;
    eqm_cf.X_s = X_s;
    eqm_cf.X = sum(X00+X01+X10+X11);
    %eqm_cf.P_sF = P_sF;
    %eqm_cf.mbar = mbar;  
    %eqm_cf.f = f;
    
    eqm_cf.ImportProbQ0 = ImportProbQ0;
    eqm_cf.ImportProbQ1 = ImportProbQ1;
    eqm_cf.EPiQ0 = EPiQ0;
    eqm_cf.EPiQ1 = EPiQ1;
    eqm_cf.ExportCutoffQ = ExportCutoffQ;
    eqm_cf.ExportProbQ = ExportProbQ;

    eqm_cf.iter_inner = iter;
    eqm_cf.err_inner = err;
    
    F = eqm_cf;

    
end